A  numerical  method  for  two-point  boundary  value  problems  with  constant  coefficients  is  devel¬ 
oped  which  is  based  on  integral  equations  and  the  spectral  integration  matrix  for  Chebyshev  nodes. 
The  method  is  stable,  achieves  superalgebraic  convergence,  and  requires  0(N  logN)  operations, 
where  N  is  the  number  of  nodes  in  the  discretization.  Although  stable  spectral  methods  have  been 
constructed  in  the  past,  they  have  generally  been  based  on  reformulating  the  recurrence  relations 
obtained  through  spectral  differentiation  in  an  attempt  to  avoid  the  ill-conditioning  introduced  by 
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1  Introduction 


r  Spectral  methods  are  now  a  popular  tool  for  the  solution  of  many  types  of  partial  differential 

[  equations.  In  the  usual  formulation,  the  basic  idea  is  to  represent  the  solution  /  by  means  of  a 

I  (truncated)  series  expansion,  and  to  compute  spatial  derivatives  of  /  by  analytic  differentiation  of 

the  series.  The  linear  map  Vs  which  takes  a  vector  of  N  function  values  {/(x;)}  to  a  vector  of 
y  derivative  values  {/'(x,)}  is  known  as  the  spectral  differentiation  matrix.  The  precise  form  of 
Vs  depends  on  the  location  of  the  points  {x<}  and  the  choice  of  the  approximating  series  For 
periodic  functions,  Fourier  series  are  used  with  the  function  tabulated  at  equispaced  nodes.  For 
i  bounded  domains,  Chebyshev  or  Legendre  series  are  used  with  the  function  tabulated  at  Chebyshev 

or  Legendre  nodes,  respectively. 

Although  spectral  differentiation  is  remarkably  accurate  in  exact  arithmetic,  there  are  a  num¬ 
ber  of  difficulties  associated  with  its  use.  Di-conditioning  of  the  matrix  with  increasing  N  fre¬ 
quently  causes  degradation  of  the  observed  precision.  Furthermore,  as  recently  demonstrated  by 
►  Trefethen  and  Trummer  for  certain  hyperbolic  problems  [6],  the  time  step  restrictions  due  to  this 

ill-conditioning  can  be  more  severe  than  those  predicted  by  the  standard  stability  theory. 

In  this  paper,  we  will  consider  only  the  simplest  steady-state  case,  namely  linear  two-point 
boundary  value  problems  with  constant  coefficients.  It  is  well  known  that  problems  of  this  type 
are  efficiently  and  accurately  solved  by  spectral  methods.  On  the  other  hand,  it  is  also  well  known 
r  that  care  must  be  taken  in  applying  spectral  methods  to  such  problems.  The  naive  approach  leads 

?  to  the  use  of  unstable  recurrence  relations  for  the  determination  of  the  expansions  coefficients,  a 

consequence  of  the  high  condition  number  of  Vs- 

|  Our  main  purpose  in  this  paper  is  to  provide  a  consistent  framework  for  developing  well- 

conditioned  spectral  methods.  After  collecting  the  necessary  results  from  approximation  theory, 
[  we  present  a  fast  algorithm  for  computing  the  indefinite  integral  of  a  given  function  by  means  of 

b  the  spectral  integration  matrix.  The  indefinite  integral  is  then  used  to  recast  the  governing  differen¬ 

tial  equation  of  the  boundary  value  problem  as  an  integral  equation,  which  is  solved  with  spectral 
accuracy. 

2  Chebyshev  Approximation 

We  will  require  several  results  from  approximation  theory.  The  Chebyshev  polynomial  of  degree  k 
on  [—1, 1]  is  defined  by  the  formula 

Tk(cos&)  =  cos (k0)  .  (1) 

.  Clearly,  |T*(x)|  <  1  for  x  €  [-1, 1], 

k 

r  T0(x)  =  1,  7\(x)  =  x  (2) 

and,  using  elementary  trigomometric  identities, 

Tk+i(x)  =  2xTk(x) -Tk-i(x)  for  k  >  1  .  (3) 

r 

The  functions  Tk  constitute  an  orthonormal  basis  with  respect  to  the  inner  product 

U,0)~  [  f(x)3(r)(  1-  x2)_1/24'x  .  (4) 

J- 1 
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The  Chebyshev  nodes  t,  of  degree  k  are  the  zeros  of  Tk ,  namely 

ti  =  cos  — ■  ,  for  i  =  0, 1,2,  ...,k  —  1  .  (5) 

Zn> 

Let  Cn[— 1,1]  denote  the  set  of  functions  defined  on  [—1,1]  with  n  continuous  derivatives.  If 
/  G  Cn[-l,l]  and 

OO 

$(*)  =  akTk(x)  (6) 

k=0 

is  the  Chebyshev  expansion  associated  with  /,  then 

ak  =  — -  f  f(x)Tk(x)(l  —  x 2)~1/2dx  =  f  /(cos  9)  cos  k6  d9  ,  (7) 

TTCfc  1  TCI fc  /o 

where  c0  =  2  and  c*  =  1  for  k  >  0.  Moreover,  the  remainder  in  truncating  the  series  at  N  term?  is 
of  the  order 

°(ivK)  “  <s) 

In  particular,  if  /  is  infinitely  differentiable,  then  the  remainder  goes  to  zero  superalgebraically 
(faster  than  any  finite  power  of  1  /N).  For  a  more  complete  discussion,  see  Gottlieb  and  Orszag  [3]. 

Remark  2.1:  In  practice,  the  Chebyshev  series  (6)  is  truncated  at  some  finite  number  of  terms, 
say  N.  By  relation  (7),  the  coefficients  ak  are  the  coefficients  of  the  Fourier  cosine  series  of  F(6)  = 
/(cos  O').  Thus,  if  /  is  tabulated  at  equispaced  points  in  0 ,  a  condition  satisfied  by  the  Chebyshev 
nodes  {t^},  we  can  obtain  all  N  coefficients  ak  by  means  of  the  FFT  or,  more  precisely,  the  Fast 
Cosine  Transform,  using  0(N  log  N)  operations.  Similarly,  the  inverse  cosine  transform  can  be  used 
to  compute  function  values  g(x)  ss  /(x)  at  the  nodes  {/, }  from  the  coefficients  ak  of  the  expansion 

2.1  Differentiation  and  Integration  of  Chebyshev  Expansions 
Definition  2.1  Let  X  be  the  space  consisting  of  infinite  sequences  of  real  numbers 

X  =  (xo,xux2,...)  ■ 

For  any  a  G  X ,  we  will  denote  by  V{a)  the  sequence  b  given  by  the  formula 

1  00 

bk  =  —  J2  P'°p  (9) 

p-fc+i 
p-f  *  odd 

where  Co  =  2  and  Ck  =  1  for  k  >  0.  We  will  refer  to  the  mapping  V  :  X  — ►  X  as  spectral 
differentiation.  W*  will  denote  by  1(a)  the  sequence  d  given  by  the  formulae 

dk  =  ^(ck-i  ■  ak-i  ~  °k+i  ■  ak+i)  for  k  >  l  (10) 

dn  =  2di  —  2d7  +  2d:> - .  (11) 

where  the  coefficients  c k  are  defined  as  above.  We  will  refer  to  the  mapping  1 :  X  — *  A"  as  spectral 
integration. 
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The  above  definitions  are  motivated  by  the  following  lemma,  which  may  be  found  in  the  Ap¬ 
pendix  to  [3], 

Lemma  2.1  Let  f  be  a  smooth  function  given  by  a  Chebyshev  series 

/(*)  =  £>*r*(x).  (12) 

*=o 

Then  the  derivative  of  f  has  a  series  expansion  of  the  form 

/'(*)  =  f>*rfc(x)  (is) 

*=o 

with  bk  given  by  (9).  The  integral  of  f  has  a  series  expansion  of  the  form 

/  f(t)dt  =J2dkTk(x)  (14) 

k= 0 

where  dk  is  given  by  (10)  and  (11). 

Remark  2.2:  The  series  expansion  (14)  is  the  basis  for  Clenshaw- Curtis  quadrature  [2].  After 
obtaining  the  coefficients  dk,  one  computes 


r  l  ~ 

/  f(x)dx  =  2  d2fc+i  . 

J- 1 


which  follows  immediately  from  the  equalities 

Tk(-l)  =  (-1)*  and  Tk(  1)  =  1  .  (16) 

It  is  clear  from  (9)  and  (10)  that  V  is  unbounded  in  the  /<»  norm,  while  HZHoo  <  2.  This 
behavior  is  reflected  in  the  condition  numbers  of  the  finite-dimensional  analogs  of  these  operators. 
If  /  is  represented  by  a  truncated  Chebyshev  expansion 

N 

f(x)=^akTk(x),  (17) 

k= o 

then  the  coefficients  of  f  are  still  given  by  (9),  but  the  summation  is  truncated  at  N  terms.  Now 
let  a  -  (a0,  a\, ... ,a jv),  a  =  (a0  +  c,  a]  -fi  e,  ...,aN  +  f),  P(a)  =  b  and  P(a)  =  b.  Then 


the  maximum  error  being  incurred  for  the  calculation  of  bo-  In  other  words,  the  process  of  differ¬ 
entiation  via  Chebyshev  series  has  a  condition  number  proportional  to  N2.  On  the  other  hand,  it 
is  easy  to  show  that  for  any  a, 

H**  ~dlj^  <  _  ,  ns) 

j|a  -  ajjoo  “  2  ’ 

where  d  =  Z(a)  and  d  =  Z(a).  In  other  words,  the  process  of  differentiation  via  Chebyshev  series 
has  a  condition  number  bounded  by  2. 
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2.2  The  Spectral  Integration  Matrix 

The  spectral  differentiation  matrix  for  Chebyshev  nodes  can  be  expressed  in  terms  of  V  by  the 
formula 

VN=Cf?-V-CN  (20) 

where  C\  is  the  discrete  cosine  transform  of  dimension  N . 

Definition  2.2  The  spectral  integration  matrix  for  Chebyshev  nodes  1/j-  is  defined  by  the  formula 

In  =  Cjv*  ‘  %  '  Cn  •  (21) 

It  is  clear  that  the  matrix  X.v  can  be  applied  to  a  vector  in  0(N  log  IV)  operations  by  using  a 
Fast  Cosine  Transform  algorithm. 

Before  turning  to  the  solution  of  two-point  boundary  value  problems,  we  briefly  investigate 
the  behavior  of  Vj<j-  and  Tn  with  a  set  of  three  examples  (Fig.  1).  Note  that  in  these  examples, 
we  test  XV  hy  differentiating  a  function  f{x)  and  we  test  Tv  by  integrating  the  corresponding 
derivative  /'(x).  The  integration  of  f(x)  itself  is  even  more  stable  and  of  less  interest.  When 
f(x)  =  sin(x),  we  observe  the  expected  convergence  as  soon  as  the  number  of  sampling  points  is 
sufficient,  approximately  tt  per  wavelength.  However,  as  the  number  of  points  increases,  differen¬ 
tiation  becomes  less  and  less  accurate  while  integration  is  essentially  unaffected.  Although  there 
is  no  need  to  use  1000  points  to  resolve  sin(x)  alone,  the  error  introduced  by  its  differentiation  re¬ 
mains  when  the  problem  becomes  more  complex  and  more  points  are  required.  For  example,  when 
/(x)  =  sin(x)  +  .Olsin(lOx),  about  50  points  are  required  to  achieve  spectral  accuracy,  at  which 
point  an  error  of  the  order  10~4  has  already  been  incurred.  When  /(x)  =  sin(x)  +  .005  sin(60x), 
the  situation  is  worse.  In  single  precision,  even  with  the  optimal  choice  of  N,  the  mean  square  error 
is  greater  than  1%. 

3  Two-point  Boundary  Value  Problems 

The  two-point  boundary  value  problems  considered  here  are  second  order  equations  of  the  form 

L  u  =  u"  +  p.\L  +  u  u  =  /(x)  ,  x  £  [—1,  l]  (22) 

with  p,  v  €  R  and  Dirichlet  conditions 


u(-l)  =  q  ,  u(l)  =  /?  .  (23) 

The  fact  that  high  order  polynomial  approximations  achieve  superalgebraic  convergence  for  such 
differential  equations  has  been  known  for  a  long  time.  Ciarlet,  Schultz  and  Varga  [1]  have  shown  that 
superalgebraic  convergence  can  be  achieved  even  when  the  governing  ordinary  differential  equation 
is  nonlinear,  so  long  as  the  solution  is  sufficiently  smooth.  In  the  standard  spectral  formulation  [3], 
using  a  “Chebyshev-tau”  method,  we  seek  a  solution 

N 

UN{x)  =  'jTakTk(x)  ,  (24) 

k= o 
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Figure  1 

The  numerical  behavior  of  spectral  differentiation  and  integration  is  demonstrated  with  three  examples. 
The  left  hand  of  each  pair  of  figures  is  a  plot  of  f(x)  as  a  dashed  curve  and  f\x )  as  a  dotted  curve.  The 
right  hand  of  each  pair  is  a  plot  of  the  mean  square  error  in  the  spectral  differentiation  of  f(x)  (dotted 
line)  and  the  spectral  integration  of  f'(x)  (dashed  line)  vs.  the  number  of  points  in  the  discretization.  For 
the  top  figure  f(x)  =  sin(x),  for  the  middle  figure,  f(x)  =  shl(x)  +  .01  sin(lOx)  and  for  the  bottom 
figure  f(x)  =  sin(x)+  .005  sin(60x). 


subject  to  the  boundary  conditions 


^2(-l)kak  =  a  and  ak  =  P  ,  (25) 

fc=o  fe=o 

where  we  have  used  the  equalities  in  (15).  From  the  differential  equation  (21)  and  the  spectral 
differentiation  matrix,  we  obtain  the  relations 

.  JV  JV 

—  V*  p(p2  -  fc2)ap  +  —  5Z  PaP  +  l/ak  =  fk  for  k  =  0,...,N-2  (26) 


P**  +  2 
p-ffc  even 


Cfc  p«k+l 
p-f  k  odd 


where  the  {/*,}  are  the  Chebyshev  expansion  coefficients  for  the  right-hand  side  f(x).  This  set  of 
equations  is  inherently  ill-conditioned.  A  reformulation  of  the  recurrence  relations  is,  therefore, 
used  to  compute  the  solution  in  a  stable  manner  (see  [3],  p.  118.). 

Consider  now  the  one-dimensional  Poisson  equation 

u"  =  f(x)  ,  ti(-l)  =  o  ,  u(l)  =  (3  .  (27) 

Rather  than  setting  up  the  recurrence  relations  as  in  (25),  it  would  be  very  attractive  to  be  able 
to  write 


u(x)  =  £  £  /(t)  dr  dt  +  CiX  -)-  Co  . 


The  spectral  integration  matrix  of  Definition  2  1  provides  us  with  precisely  the  ability  to  compute 
this  formal  solution.  The  constants  C\  and  Co  are  chosen  to  satisfy  the  boundary  condition. 

When  the  differential  equation  contains  terms  in  u'  or  u,  we  still  have  a  simple  analytic  expres¬ 
sion  for  the  solution.  We  may  assume,  without  loss  of  generality,  that  we  are  given  homogeneous 
boundary  conditions  and  that  the  corresponding  Green’s  function  has  the  form 


Uj(x)  v\(t)  for  x  <  t 
U2(x)v2(t)  for  x  >  t 


where  u}  and  U2  are  solutions  of  the  homogeneous  equation  Lu  —  0.  Since  the  equation  has 
constant  coefficients,  Ui  and  U2  are  known  explicitely.  The  desired  solution  can  then  be  written  as 

u(x)  =  J^G(x,t)f(t)dt  (30) 

=  Uj(x)  J  j  v\{t)f{t)dt  +  (31) 

«2(x)  •  (y^  V2(t)f(t)  dt  -  £  V2{t)f{t)  dt)  .  (32) 

The  indefinite  integrals  in  the  preceding  expression  can  be  tabulated  by  means  of  the  spectral 
integration  matrix,  while  the  definite  integral  can  be  computed  by  formula  (14).  Once  this  initial 
work,  requiring  two  Fast  Cosine  Transforms,  is  done,  the  solution  is  obtained  using  approximately 
3.V  additional  operations,  where  N  is  the  number  of  Chebyshev  nodes  used.  A  similar  observation 
is  made  by  Rokhlin  in  [4],  where  the  integrals  are  computed  by  a  finite-order  quadrature  formula. 
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His  approach  has  the  advantage  that  it  does  not  denend  on  the  location  of  the  discretization  nodes. 
One  drawback  of  this  particular  use  of  Green’s  functions  is  that  the  terms  uy,  t»lt  U2  and  v2  may 
behave  much  more  violently  than  either  the  right-hand  side  or  the  solution,  requiring  many  more 
points  in  the  discretization  than  necessary.  Another  drawback  is  that  the  method  does  not  extend 
to  non-constant  coefficient  problems,  since  we  cannot  determine  the  Green’s  function  analytically. 

We  consider  a  closely  related  approach  which  does  not  require  knowledge  of  the  Green’s  function. 
An  integral  equation  is  constructed  by  solving  for  a(x)  =  u"(x )  rather  than  u  itself.  The  original 
system  (21)  becomes 

a(x)-4-/.tJ  <r(t)  dt  +  nCi  +  v J  J  ct(t)  dr  dt  +  v  Cyx  +  u  C0  =  f(x)  (33) 

Representing  o{x )  and  f(x)  by  truncated  Chebyshev  series 


v(x)  =  5Z  anTn(x) 


f(x)  =  '£fnTn(x) 


and  using  the  spectral  integration  matrix,  we  obtain  the  system  of  equations 

ao  +  n  C\  +  v  Co  =  fo 
ai  +  ^(co  «o  ~  C2  <*2)  4-  -(8  Ci  +  aj  -  a3)  =  /1 


f.1  1/ 

ak  +  OL^Ck~l  afc_1  ~  Cfc+1  afc+i)  +  (c*-i  dk-i  ~  ck+\  dt+i)  =  fk 


(37) 

for  k  =  2, ...,  N  ,  (38) 


where  c0  =  2,  Ci  =  ...  =  cp?  =  1,  c/t  =  0  for  k  >  N,  and  the  coefficients  dk  are  given  by  equation 
( 10).  This  is  a  system  of  N  + 1  equations  with  N  +  3  unknowns  (the  coefficients  a*  and  the  constants 
of  integration  Co  and  Cy).  Two  additional  equations  are  obtained  from  the  boundary  conditions 

N  1 

Co-Ci-f^— (c^dfc-x-Cfc+id^X-l)*  =  a,  (39) 

k—1 

N  1 

Co  +  Ci  +  7^r(ck-i  dk-i  —  ck+i  dk+i)  =  0  ■  (40) 

k= 2 

The  discrete  problem  is  pentadiagonal  except  for  the  two  rows  derived  from  the  boundary  condi¬ 
tions,  and  can  be  solved  using  approximately  ION  floating  point  operations. 


4  Numerical  Examples 

A  two-point  boundary  value  problem  solver,  using  the  method  of  the  previous  section,  has  been 
implemented  and  tested  on  a  variety  of  examples.  It  requires  two  cosine  transforms  and  the  solution 
of  one  linear  system,  with  a  total  computational  cost  estimated  at  10  A  (log  A  +  1)  floating  point 
operations.  All  calculations  cited  below  were  carried  out  in  double  precision  on  a  SUN  3/50 
workstation  with  f68881  floating  point  accelerator. 


Figure  2 

The  numerical  behavior  of  the  integral  equation  algorithm,  using  three  examples.  The  left  hand  of  each 
pair  of  figures  is  a  plot  of  the  exact  solution  as  a  function  of  X.  The  right  hand  of  each  pair  is  a  plot  of  the 
mean  square  error  in  the  computed  solution  vs.  the  number  of  points  in  the  discretization.  The  governing 
differential  equations  are  discussed  in  the  text. 
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•  Table  1 

Table  of  CPU  times  (in  seconds)  and  mean  square  errors  in  computing  the  solutions  to  three  boundary 
value  problems.  The  subscripts  1-3  refer  to  the  three  equations  discussed  in  the  text  and  displayed 
graphically  in  Fig.  2  from  top  to  bottom. 
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4 


4 
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The  behavior  of  the  algorithm  is  demonstrated  with  three  examples  (Fig.  2  and  Table  1  ).  In 
the  first  case,  we  used  a  model  problem  from  Stoer  and  Bulirsch  [5] 

-  y"  +  400y  =  -400  cos2  xx  —  27T2  cos  2xx  ,  (41) 


y(0)  =  y(i)  =  o  , 


with  exact  solution 


y(x)  - 


e20x  + 


e  201  -  cos2  xx 


Standard  finite  difference  and  finite  element  methods  tend  to  converge  quite  slowly,  due  to  the 
large  derivatives  of  the  exact  solution  near  the  boundaries.  While  multiple  shooting,  which  is  rec¬ 
ommended  in  [5],  is  a  viable  approach,  the  method  is  computationally  expensive.  Our  calculations 
show  that  the  corresponding  integral  equation  is  solved  to  spectral  accuracy  with  very  little  efTort. 

The  second  example  involves  a  boundary  layer  near  each  endpoint.  The  governing  equation  is 

cy"-y=  0,  (44) 


y(— i)  =  i ,  y(i)  =  2,  (45) 

where  e  =  10“ J.  As  is  well-known,  the  Chebyshev  nodes  are  particularly  good  at  resolving  boundary 
layers  since  they  tend  to  cluster  at  the  two  endpoints  (see  [3]). 

The  third  example  is  one  where  the  solution  is  very  oscillatory: 

y"  +  5  y'  +  10000  y  =  -500  cos( lOOx)  e"5x  ,  (46) 

y(0)  =  0  y(l)  =  sin(  100)  e~5  ,  (47) 

for  which  the  exact  solution  is 

y(x)  =  sin(100x)  e'5x  .  (48) 

In  each  case,  it  is  clear  that  the  spectral  integral  formulation  is  both  rapidly  convergent  and 
stable. 
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5  Conclusions 


The  spectral  integration  matrix  is  a  well-conditioned  operator  which  yields  an  antiderivative  of 
a  function  tabulated  at  Chebyshev  nodes.  In  this  paper,  we  have  presented  a  fast  algorithm  for 
the  solution  of  constant  coefficient  two-point  boundary  value  problems  through  the  use  of  integral 
equations  and  spectral  integration.  The  difficulty  with  variable  coefficient  problems  lies  not  in  the 
formulation  of  the  integral  equation,  but  in  the  fact  that  the  resulting  system  of  equations  for 
the  coefficients  of  the  Chebyshev  series  of  n"  is  dense.  Gaussian  elimination  would  require  0(.V3) 
operations,  where  N  is  the  number  of  Chebyshev  nodes  used  in  the  discretization.  On  the  other 
hand,  the  spectral  integration  matrix  can  be  used  to  apply  the  integral  operator  in  0(N  log  N) 
operations,  making  iterative  methods  much  more  attractive.  Rokhlin  [4]  haa  demonstrated  that 
conjugate  residual  type  methods  can  work  quite  well  for  such  integral  equations.  The  number  of 
iterations  required  is  a  function  of  the  underlying  problem,  and  does  not  increase  with  the  number 
of  nodes.  Unfortunately,  for  many  situations  of  interest,  complex  behavior  of  the  solution  causes 
the  condition  number  of  the  underlying  problem  and  the  number  of  iterations  to  be  large,  so  that 
direct  methods  would  be  preferable. 


I  would  like  to  thank  V.  Rokhlin  and  F.  Saied  for  many  useful  discussions. 
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